!     Copyright (C) 2010 Benjamin Piaud
!
!     LIMBES is free software; you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation; either version 3, or (at your option)
!     any later version.

!     LIMBES is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with LIMBES; if not, see <http://www.gnu.org/licenses/>


program LIMBES_mesh
  use LIMBES_mesh_mod
  implicit none
  integer*4 :: ni,n1,n2,xi,yi,x1,y1,wall
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!variables can be added by user!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
  integer*4 :: cr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!domain size and delta x indicate by user!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  nx=1002                ! numbers of cells in x axis
  ny=300              ! numbers of cells in y axis
  ! for output VTK files
  VTK_xmin=0         ! xmin of interesting domain to write output result
  VTK_xmax=0         ! xmax of interesting domain to write output result
  VTK_ymin=0         ! ymin of interesting domain to write output result
  VTK_ymax=0         ! ymax of interesting domain to write output result
  delta_x=5.0d-8     ! delta x
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


  Nmax=nx*ny
  allocate(phase(1:nx,1:ny))  
  allocate(cell_xi(1:Nmax))
  allocate(cell_yi(1:Nmax))
  allocate(cell_ni(1:nx,1:ny))


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!mesh defined by user!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! user must to indicate if the cell at position xi,yi is included in the fluid domain
  ! by writing phase(xi,yi)=1
  ! if phase(xi,yi)=0 that means, this cell is in the "solid"
  ! if the cell at xi=1 or xi=nx or yi=1 or yi=ny is the fluid domain so the boundary
  ! condition will be periodic. 
  do xi=1,nx
     do yi=1,ny
        phase(xi,yi)=1
        ! left side and right side
        if (xi.EQ.1.OR.xi.EQ.nx) phase(xi,yi)=0
        ! zone 1
        if (xi.GT.1.AND.xi.LE.495.AND.yi.LE.100) phase(xi,yi)=0
        ! zone 2
        if (xi.GE.506.AND.xi.LT.nx.AND.yi.LE.100) phase(xi,yi)=0
        ! zone 3
        if (xi.GT.1.AND.xi.LE.495.AND.yi.GE.201) phase(xi,yi)=0
        ! zone 4
        if (xi.GE.506.AND.xi.LT.nx.AND.yi.GE.201) phase(xi,yi)=0
     enddo
  enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


  ni=0
  do xi=1,nx
     do yi=1,ny
        if  (phase(xi,yi).EQ.1) then
           ni=ni+1
           cell_xi(ni)=xi
           cell_yi(ni)=yi
           cell_ni(xi,yi)=ni
        endif
     enddo
  enddo
  fluid_ntot=ni

  ! detecting boundary 
  call boundary()
  ! determine cell neighbours
  call neighbour()

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!boundary properties defined by user!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! for each wall, you must indicate boundary type (boundary_type) and 4 properties (alpha(1or2)) 
  ! the boundary type can be:
  ! 'S' for specular reflexion. No properties are needed you can write alpha(ni,1 to 4)=zero.
  ! 'T' for diffuse reflexion (or Thermal condition).
  !      alpha(ni,1) is accomodation coefficient (between 0.0 and 1.0) 
  !      alpha(ni,2) is the tangential wall velocity
  !      alpha(ni,3) is the wall temperature
  !      alpha(ni,4)=zero (not used)
  ! 'F' for imposed flow (to represent open system with the Flow known at the infinite)
  !      alpha(ni,1) is the density of the flow
  !      alpha(ni,2) is the x component of the velocity
  !      alpha(ni,3) is the y component of the velocity
  !      alpha(ni,4) is the temperature
  ! 'PT' for imposed Pressure and Temperature (to represent a controled pool)  
  !      alpha(ni,1) is the pressure of the pool
  !      alpha(ni,2) is the temperature of the pool
  !      alpha(ni,3)=zero (not used)
  !      alpha(ni,4)=zero (not used)
  ! 'P' for imposed Pressure (to represent a controled pool)  
  !      alpha(ni,1) is the pressure of the pool
  !      alpha(ni,2)=zero (not used)
  !      alpha(ni,3)=zero (not used)
  !      alpha(ni,4)=zero (not used)
  do ni=1,fluid_ntot
     xi=cell_xi(ni)
     yi=cell_yi(ni)
     wall=wall_RLUD(ni,1)+ wall_RLUD(ni,2)+ wall_RLUD(ni,3)+ wall_RLUD(ni,4)
     if (wall.EQ.1) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!boundary properties for cell with 1 wall!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! wall of the main channel
        if (yi.EQ.101.OR.yi.EQ.200) then
           boundary_type1(ni)='T'
           alpha1(ni,1)=1.0D+0
           alpha1(ni,2)=0.0D+0
           alpha1(ni,3)=300.0D+0
           alpha1(ni,4)=zero
        endif
        ! wall of the second channel
        if (xi.EQ.496.OR.xi.EQ.505) then
           boundary_type1(ni)='T'
           alpha1(ni,1)=1.0D+0
           alpha1(ni,2)=0.0D+0
           alpha1(ni,3)=330.0D+0
           alpha1(ni,4)=zero
        endif
        !left side
        if (xi.EQ.2) then
           boundary_type1(ni)='PT'
           alpha1(ni,1)=1.179D+0*287.0D+0*300.0D+0
           alpha1(ni,2)=300.0D+0
           alpha1(ni,3)=zero
           alpha1(ni,4)=zero
        endif
        !left side
        if (xi.EQ.nx-1) then
           boundary_type1(ni)='PT'
           alpha1(ni,1)=1.177D+0*287.0D+0*300.0D+0
           alpha1(ni,2)=300.0D+0
           alpha1(ni,3)=zero
           alpha1(ni,4)=zero
        endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!end of boundary properties for cell with 1 wall!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     endif
     if (wall.EQ.2) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!boundary properties for cell with 2 walls!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !left side
        if (xi.EQ.2) then
           boundary_type2(ni,1)='PT'
           alpha1(ni,1)=1.179D+0*287.0D+0*300.0D+0
           alpha1(ni,2)=300.0D+0
           alpha1(ni,3)=zero
           alpha1(ni,4)=zero
           boundary_type2(ni,2)='T'
           alpha2(ni,1)=1.0D+0
           alpha2(ni,2)=0.0D+0
           alpha2(ni,3)=300.0D+0
           alpha2(ni,4)=zero
        endif
        !right side
        if (xi.EQ.(nx-1)) then
           boundary_type2(ni,1)='PT'
           alpha1(ni,1)=1.177D+0*287.0D+0*300.0D+0
           alpha1(ni,2)=300.0D+0
           alpha1(ni,3)=zero
           alpha1(ni,4)=zero
           boundary_type2(ni,2)='T'
           alpha2(ni,1)=1.0D+0
           alpha2(ni,2)=0.0D+0
           alpha2(ni,3)=300.0D+0
           alpha2(ni,4)=zero
        endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!end of boundary properties for cell with 2 walls!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     endif
  enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !writing output
  call write_mesh()
  call write_boundary()

end program LIMBES_mesh
